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Abstract. This paper reports on a significant advance in the area of non- 
reflecting boundary conditions (NRBCs) for unsteady flow computations. As 
a part of the development of the space-time conservation element and solution 
element (CE/SE) method, sets of NRBCs for ID Euler problems are developed 
without using any characteristics- based techniques. These conditions are much 
simpler than those commonly reported in the literature, yet so robust that 
they are applicable to subsonic, transonic and supersonic flows even in the 
presence of discontinuities. In addition, the straightforward multidimensional 
extensions of the presnt ID NRBCs have been shown numerically to be equally 
simple and robust. The paper details the theoretical underpinning of these 
NRBCs, and explains their unique robustness and accuracy in terms of the 
conservation of space-time fluxes. Some numerical results for an extended 
Sod’s shock-tube problem, illustrating the effectiveness of the present NRBCs 
are included, together with an associated simple Fortran computer program. 
As a preliminary to the present development, a review of the basic CE/SE 
schemes is also included. 


1. Introduction 

Because of computing resource limitation and other considerations, it is often 
required that the computational domain used in a flow simulation represent only 
a subdomain of a larger physical domain. To obtain a numerical solution that 
closely resembles the physical flow field in this subdomain, ideally the conditions at 
the computational boundary should be specified using the physical flow conditions 
there. Unfortunately, these conditions generally are not known without first solving 
the larger physical flow field. 

Despite the above difficulty, with proper boundary treatments, an accurate 
simulation of the physical flow over a subdomain using a computational domain 
that covers only the subdomain is possible if certain conditions are met. As an 
example, assume that flow disturbances are generated within the subdomain while 
no disturbances enter it from outside. For this case, an accurate numerical solution 
over the subdomain can be obtained by imposing proper nonreflecting conditions at 
the computational boundary. These nonreflecting boundary conditions (NRBCs) 
are designed such that flow disturbances can propagate out of the computational 


NASA/TM— 2003-2 1 2495/REV 1 


1 


domain smoothly without inducing substantial spurious reflections from the bound- 
ary. Such reflections can distort the computed solution and render it completely 
useless. 

Design and application of NRBCs have been a research topic for a long time [1- 
8]. In the following, we describe three main established approaches for implemeting 
NRBCs. 

The first approach is based on ID characteristic decomposition of flow vari- 
ables and it is most suitable w r hen the waves propagate toward the boundary in 
the normal direction. Engquist and Majda [1] express the nonreflecting boundary 
condition as the requirement that the local perturbation propagating along the 
incoming characteristics be made to vanish. In practice, this consists of project- 
ing the flow equations onto the normal direction of the boundary, converting the 
conservative variables to characteristic variables, finding the characteristics that 
enter the domain, and finally setting the corresponding characteristic variables to 
zero. Although such a procedure is delicate and tedious, it still represents the most 
commonly used nonreflecting boundary treatment [2,3,5] 

The second approach is based on the far-field asymptotic solutions [4] and it is 
ideal when the mean flow near the boundary is nearly uniform. 

In the third approach [6 8], efforts are made to insure that the disturbances in 
the buffer zone that lies outside the computatinal boundary will not reflect back 
into the computational domain. The most recent development in this area is the 
so called perfectly matched layer (PML) method [7]. 

Generally speaking, the established methods described above are not designed 
for flow problems involving shocks and contact discontinuities. Furthermore, their 
implementation generally requires a considerable effort. 

As an integral part of the development of the space-time conservation element 
and solution element (CE/SE) method [9 40], several sets of NRBCs for ID Euler 
problems will be derived in this paper using a nontraditional concept based entirely 
on a simple assumption about the space-time flux distribution in the neighborhood 
of a spatial boundary. As such, the derivation of these NRBCs is consistent with 
the two basic tenets [15, p. 89] of the CE/SE development, i.e., (i) to capture 
physics more efficiently and realistically, the modeling focus should be placed on 
the original integral form of the physical conservation laws, rather than the differ- 
ential form (which follows from the integral form under the additional smoothness 
assumption); and (ii) to simplfy mathematics, the use of any elaborate partial dif- 
ferential equation theory should be avoided as much as possible. As will be shown, 
the derived NRBCs are indeed very simple (e.g., the solution values at a boundary 
mesh point are simply taken from those at a neighboring interior mesh point), and 
yet so robust that they are applicable to subsonic, transonic and supersonic flows 
even in the presence of discontinuities. In addition, as will be described further in 
Sec. 6, the straightforward multidimensional extensions of the present ID NRBCs 
have been shown numerically to be equally effective and robust. 

The space-time CE/SE method is a high-resolution and genuinely multidimen- 
sional method for solving conservation laws. It is not an incremental improvement 
of a previously existing method, and it has many nontraditional features. They 
include: (i) a unified treatment of space and time, (ii) the introduction of con- 
servation elements (CEs) and solution element (SEs) as the vehicles for enforcing 
space-time flux conservation, and (iii) a time marching strategy that has a space- 
time staggered stencil at its core and, as such, can capture shocks without using 
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Riemann solvers. Note that conservation elements are nonoverlapping space- time 
subdomains introduced such that (i) the computational domain can be filled by 
these subdomains; and (ii) fiux conservation can be enforced over each of them 
and also over the union of any combination of them. On the other hand, solution 
elements are nonoverlapping space-time subdomains introduced such that (i) the 
boundary of any CE is covered by a combination of SEs; and (ii) any physical flux 
vector is approximated using simple smooth functions within a SE. In general, a 
CE does not coincide with a SE. 

Without using preconditioning or other special techniques, since its inception 
[9] the CE/SE method has been used to obtain numerous highly accurate ID, 
2D and 3D steady and unsteady flow solutions with Mach numbers ranging from 
0.0028 to 10. The flow phenomena modeled include traveling and interacting shocks, 
acoustic waves, shedding vortices, shock/boundary-layer interaction, detonation 
waves, cavitation and hydraulic jump. In particular, the rather unique capability 
of the CE/SE method to resolve both strong shocks and small disturbances (e.g., 
acoustic waves) simultaneously [19,21 26] makes it a unique tool for attacking the 
problems in computational aeroacustics (CAA). Note that the fact that the (second 
order) CE/SE method can solve CAA problems accurately is an exception to the 
commonly-held wisdom that a second-order scheme is not adequate for solving 
CAA problems. Also note that, while numerical dissipation is required for shock 
capturing, it may also result in annihilation of small disturbances. Thus a solver 
that can handle both strong shocks and small disturbances simultaneously must be 
able to overcome this difficulty. 

The rest of the paper is organized as follows: A brief description of the CE/SE 
method is provided in Sec. 2. The concept of generalized conservation elements 
is introduced in Sec. 3. Using the preliminaries given in Secs. 2 and 3, several 
sets of NRBCs are derived in Sec. 4. Numerical results that validate the derived 
NRBCs are given in Sec. 5. The concluding remarks are given in Sec. 6. Finally, 
to give the reader a concrete example about the simplicity and robustness of the 
CE/SE method in general and the the current NRBCs in particular, a short Fortran 
program for solving an extended Sod’s shock tube problem with a set of NRBCs 
imposed at its two open ends is listed in Appendix. 


2. Review of the CE/SE method 


As an example, the ID CE/SE method will be described by considering the 
PDE 


( 2 . 1 ) 


du du 

m +a d^-° 


where a is a constant. Let x\ = x, and x 2 = t be considered as the coordinates of 
a two-dimensional Euclidean space E>. Then, by using Gauss’ divergence theorem 
in the space-time E 2 , it can be shown that Eq. (2.1) is the differential form of the 
integral conservation law 


( 2 . 2 ) 


(p h • ds = 0 
JS(V) 


As depicted in Fig. 1, here (i) S(V) is the boundary of an arbitrary space-time 
region V in E 2 , (ii) h = (au, u), and (iii) ds = dan with da and n, respectively, 
being the area and the unit outward normal of a surface element on S(U). Note 
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that: (i) because h • ds is the space- time flux of h leaving the region V through 
the surface element ds. Eq. (2.2) simply states that the total space-time flux of h 
leaving V through S(V) vanishes; (ii) in E 2 , da is the length of a line segment on 
the simple closed curve S(V); and (iii) all mathematical operations can be carried 
out as though E 2 were an ordinary two-dimensional Euclidean space. 

To proceed, let ft denote the set of all space-time staggered mesh points (j,n) 
in E 2 (dots in Fig. 2(a)), where n = 0 , ±1/2, ±1, ±3/2, ±2, . . ., and, for each n, 

j = n ± 1/2, n ± 3/2, n ± 5/2, Each ( j,n ) G 12 is associated with a solution 

element, i.e., SE O’, n). By definition, SE(j, n) is the interior of the space-time 
region bounded by a dashed curve depicted iri Fig. 2(b). It includes a horizontal 
line segment, a vertical line segment, and their immediate neighborhood. 

For any (x,t) G SE O’, n), u(x,t) and h(x,t), respectively, are approximated by 

(2.3) u*(x, t ,j, n) = f u? + (u x )](x-Xj) + (u t )](t - t n ) 
and 

(2.4) h m (x,t;j,n) d = ( au m (x,t;j,n ), u’(x,t;j,Ji)) 


Note that (i) u ", ( u x ) T j , and ( u t ) T j are constants in SE(j, n), (ii) (xj , t n ) are the 
coordinates of the mesh point (j, n), and (iii) Eq. (2.4) is the numerical analogue 
of the definition h = (au, u). 

Let u = u*(x,t;j,n) satisfy Eq. (2.1) within SE(j, n). Then one has (u*)t l = 
—a(u x )y As a result, Eq. (2.3) reduces to 

(2.5) u*(x,t;j,n)=u" + {u x )j[(x-x j )-a(t-t n )], (x,t) G SE(j,n) 

i.e., u n j and ( u x )j are the only independent marching variables associated with 

U,n). 

Let E 2 be divided into nonoverlapping rectangular regions (see Fig. 2(a)). As 
depicted in Figs. 2(c)-2(e), (i) two such regions, i.e., CE -(i, n) and CE-j.^', n), are 
associated with each interior mesh point (j,n) G f2; and (ii) CE(j,n) is the union 
of CE -(j,n) and CE +(j,n). 

Note that, among the line segments forming the boundary of CE_ (j, n), AB and 
AD belong to SE^n), while CB and CD belong to SE (j - 1/2 ,n — 1/2). Similarly, 
the boundary of CE+(j, n) belongs to either SE(j, n) or SE(j + 1/2, n — 1/2). As a 
result, by imposing two conservation conditions at each (j,n) G ft, i.e., 


(2.6) (p h* • ds — 0, and (p h* • ds = 0, (j,n) G ft 

J S(CE+(j,n)) Js(CE-(j,n )) 

and using Eqs. (2.4) and (2.5), one can obtain two equations for the two unknowns 
u ™ and {u x )j . By solving these equations, one has (i) 


(2.7) 


»? = \ {a + ^ + (i - + o - o [k + )":.« - «)«’/=] } 

and, assuming 1 — v 1 ^ 0, (ii) 


( 2 . 8 ) («+)? = - [ vk $ - - (1 - - (1 + ^)«)” +1 ‘ / / 2 


+ vn-l/2 


,+^n-l/2l 


J 
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Here, (i) v ( = f aAt/Ax, and (ii) for any (j,n) E fi, 

(2-9) («+)? d =^(«x)!? 

Note that derivation of Eqs. (2.7) and (2.8) can be facilitated by the following 
observations: because u*(x,t;j,n) is linear in x and t , it can be shown that the 
total flux of h* leaving CE_(j, n) or CE+(j, n) through any of the four line segments 
that form its boundary is equal to the scalar product of the vector h* evaluated at 
the midpoint of the line segment and the “surface” vector (i.e., the unit outward 
normal multiplied by the length) of the line segment. 

The a scheme [9,11,17], the explicit nondissipative CE/SE solver for Eq. (2.1), 
is formed by Eqs. (2.7) and (2.8). Because, for any (j, n) E H, the total flux of h* 
leaving each of CE_(j,n) and CE +(j,n) vanishes (see Eq. (2.6)), CE-(j,n) and 
CE+ (j, n), (j, n) E il, will be referred to as the conservation elements (CEs) of the 
a scheme. In addition, because the surface integration over any interface separating 
two neighboring CEs is evaluated using the information from a single SE, obviously 
the flux leaving one of these CEs through the interface is the negative of that leaving 
another CE through the same interface. As a result, the local conservation relations 
Eq. (2.6) lead to a global flux conservation relation, i.e., the total flux ofh * leaving 
the boundary of any space-time region that is the union of any combination of 
CEs will also vanish. In particular, because CE(j, n) is the union of CE_ (j, n) and 
CE+(j,n), 

(2.10) (f h*-ds=0, (j,n)€fi 

must follow from Eq. (2.6). In fact, it can be shown that Eq. (2.10) is equivalent 
toEq. (2.7). 

In addition to the nondissipative a scheme, there is a family of dissipative 
CE/SE solvers of Eq. (2.1) in which only the less stringent conservation condition 
Eq. (2.10) is assumed [11]. Because Eq. (2.10) is equivalent to Eq. (2.7), for each of 
these schemes, u ” is still evaluated using Eq. (2.7) while (w+ )” is evaluated in terms 

°fu?-% 2 and ('ut)j± i /2 us hig an equation different from Eq. (2.8). Hereafter any 
member of this family of dissipative extensions of the a scheme will be referred to as 
an a 1 scheme. Among the a' schemes is one (referred to as the a-a scheme) which is 
among the simplest and yet capable of handling solutions with discontinuities. For 
this scheme, ( u +)” is evaluated using a finite-difference/ weighted-average procedure 
which involves a parameter a (see Eqs. (2.62), (2.63) and (2.65) in [17]). Note that, 
because only Eq. (2.10), but not Eq. (2.6), is satisfied by an a 1 scheme, the CEs of 
an a' scheme are CE(j,n ), (j, n) E fi, rather than CE±(j,n), (j, n) E f l. 

The a scheme and its extensions described above have been extended to become 
solvers of systems of conservation laws in one, two and three spatial dimensions [10- 
12, 15, 1G, 20, 32]. In particular, the method was extended to two and three spatial 
dimensions by using triangles and tetrahedrons, respectively, as the basic building 
blocks of the spatial meshes [15,16,32]. In addition, it was also extended to two 
and three spatial dimensions by using quadrilaterals and hexahedrons, respectively, 
as the basic building blocks of the spatial meshes [20]. As a preliminary for Secs. 
3-5, the ID Euler equations along with the associated CE/SE solvers are briefly 
discussed in the following. 
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The ID Euler equations of a perfect gas can be expressed as [11] 

/ n i i \ dU-m dfm n 100 

(2 ' U) ■0T + &r = °’ rn= 1,2,3 

where (i) the conservative How variables u m , m = 1,2,3 are functions of the mass 
density /?, the fluid velocity v and the static pressure/?; and (ii) / m , rn = 1,2,3, are 
functions of u m , m = 1,2,3. The integral form of Eq. (2.11) in space-time E 2 is 

(2.12) (fi h m • ds = 0, m = 1,2,3 

JS{V) 

where h m = (/ m , u m ), m = 1, 2,3, are the space-time mass, momentum, and energy 
current density vectors, respectively. 

The procedures used to construct the a scheme and its extensions was extended 
to construct several ID Euler solvers [10-12]. For these Euler solvers, the indepen- 
dent unknowns at each (j,n) € Q are (ti m )j and (u mx )j, m = 1,2,3. Also, (i) the 
Euler version of Eq. (2.9) is 

(2.13) (4)? = jMJ 

(ii) the ID Euler version of the a scheme is constructed using the conservation 
conditions 

(2.14) <p h* n • ds = 0, and (p h* n • ds = 0, m = 1, 2, 3 

^S(C£ + (j,n)) Js(CE-(j,n)) 

where /i^ is the Euler version of /i*; (iii) the ID Euler version of an a' scheme 
(referred to hereafter as a ID Euler a ' scheme) is only required to satisfy the less 
stringent conservation conditions 

(2.15) <f h* n -ds = 0, m — 1,2,3 

Js(CEU.n)) 

and (iv) the ID Euler version of the a-a scheme (referred to hereafter as the ID 
Euler a-a scheme), which is defined by Eqs. (4.24) and (4.38) in [11], has also been 
shown to be an accurate and robust shock-capturing solver. 

3. Generalized Conservation Elements 

For any a' scheme, CE (j,n) (see Fig. 2(e)) is a CE. In other words, Eq. (2.10) 
is valid if, with the aid of Eqs. (2.4), (2.5) and (2.9), 

(a) h* at any point on CB and CD is evaluated using { u tYjZljo, 

(b) h* at any point on ED and EF is evaluated using an d ( w x ^+ 1/2 ’ 

and 

(c) h* at any point on AB and AF is evaluated using u 7 j and (n+)™. 

However , CE+(j, n) and CE-(j,n) (see Figs. 2(c) and 2(d)) are not CEs of an a' 
scheme. In other words, Eq. (2.6) is not valid if 

(d) h* at any point on AD (which is a part of SE O’, n)) is evaluated using 11 ™ 
and (u+)". 

Assume rules (a)-(c). Then the above discussion implies that, for an a 1 scheme, 
the total flux of h* leaving the boundary of CE_(j,n) vanishes only if we ignore 
the rule (d) and instead assume the rule: 
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(e) the flux leaving CE ~(j,n) through AD is defined to be the value such that 
the total flux leaving the boundary of CE ~(j,n) vanishes. 

Note that, unlike the rule (d), the new rule (e) does not specify h* at each individual 
point on AD. The rule (e) only defines a scalar, i.e., the total flux passing through 
the entire line segment AD. Because the flux-conservation property of CE_(j,n) 
is established by imposing the rule (e), a space-time region such as CE_ (j, n) will 
be referred to as a generalized conservation element (GCE). 

To proceed, without exception hereafter the flux leaving a space-time region A 
through an interface separating A and a neighboring region B is assumed to be the 
same flux as that entering B through the interface. As an example, the flux leaving 
CE_(j,n) through AD is the same flux that entering CE+(j,n) through AD (i.e., 
the flux leaving CE_(j, n) through AD is the negative of the flux that leaving 
CE +(j,n) through AD). Then because (i) CE (j,n) is the union of CE_(j, n) and 
CE+(j, n), and (ii) CE (j,n) is a CE for any a' scheme, the rules (a)-(c) and (e) 
lead to the conclusion that CE + (j,n) is also a GCE for any a' scheme. In fact, 
assuming the rules (a)-(c), the flux at AD defined by the rule (e) is identical to 
that defined by the rule: 

(f) the flux leaving CE +(j,n) through AD is defined to be the value such that 
the total flux leaving the boundary of CE^^', n) vanishes. 

At this juncture, note that the concept of GCEs is introduced such that the physical 
conservation law Eq. (2.2) can still have a numerical analogue in a space-time region 
that is a subset of a CE of any a' scheme. 

By using a general rule that is similar to the rule (e), we can define other GCEs 
for any a' scheme. As an example, consider Fig. 3(a). Here P denotes any point 
on BC ( P may coincide with point B or point C). Because the fluxes leaving the 
triangle ABP through BP and AB, respectively, can be evaluated using the rules 
(a) and (c), one can define the flux leaving ABP through AP to be the value such 
that the total flux leaving the boundary of ABP vanishes. With this definition, 
ABP becomes a GCE. Because the quadrilateral ABCD, (i.e., CE_(j, n)) is a 
GCE, using an argument similar to that was used to establish the equivalence of 
the above rules (e) and (f), one concludes that the quadrilateral APCD is also a 
GCE. To emphasize the fact that h* at each individual point on either AD or AP 
is not specified , AD and AP are denoted by dashed lines in Fig. 3(a). The same 
convention will also be used in the rest of this paper. 

In a similar fashion, one can divide CE_(j, n) or CE + (^, n) into other combi- 
nations of two GCEs (see Fig. 3(b)). 

Furthermore, by (i) dividing the quadrilateral APCD (a GCE) in Fig. 3(a) 
into two triangles APD and DPC (see Fig. 3(c)), and (ii) defining the flux leaving 
DPC through DP to be the value such that the total flux leaving the boundary of 
DPC vanishes, one can divide CE into three GCEs, i.e., the triangles ABP , 
APD and DPC. In a similar fashion, one can divide CE_(j,n) or CE+^n) into 
other combinations of three or more GCEs. The discussion of GCEs for any a' 
scheme is concluded with the following remarks: 

(a) Consider Fig. 4. Let P be an interior point on AD. Because the rule (e) 
does not specify the flux passing through AP or PD, one may not assign a 
flux on CP such that the triangle PCD and the quadrilateral ABCP can 
be considered as GCEs. 
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(b) Consider Fig. 5(a). Let h* be given at every point on DC , CD, DC, and GF 
(these line segments are denoted by solid lines). As a result, the flux on DC , 
CE , EA, AD, DG and GF are also given. It follows that one can assign 
a flux on AD and also a flux on EF such that both the triangle ADC and 
the quadrilateral DEFG are GCEs. Because the fluxes at all line segments 
that form the boundary of the union of ADC and DEFG are well defined 
and the total flux leaving the boundary vanishes, the union is also a GCE. 

(c) Consider Fig. 5(b). This figure is identical to Fig. 5(a) except that h* is now 
given at every point on AD, DC, DG, GF, and FE. As a result, one can 
assign a flux on AC and also a flux on DE such that both ADC and DEFG 
are GCEs. However, the question of whether the union of ADC and DEFG 
is a GCE is an ill-posed problem because (i) the flux at CE and that at AD 
are not uniquely defined, and (ii) CE and AD form part of the boundary of 
this union. 

Finally, note that the discussions of GCEs given above are applicable to any Euler 
version of any a' scheme if h* is replaced by each h* n , m = 1, 2, 3. This fact becomes 
apparent if one recalls that Eq. (2.15) is satisfied by any Euler version of any a' 
scheme, and also observes that Eq. (2.15), for each m = 1,2,3, becomes Eq. (2.10) 
if h* m is replaced by h* . 


4. Nonreflecting Boundary Conditions 

In this section, several sets of XRBCs for the CE/SE method will be derived 
using a nontraditional approach. As a prerequisite, a set of basic concepts will first 
be elaborated using a simple initial-value problem as an example. 


4.1. Basic Concepts. Consider the initial- value problem defined by Eq. (2.1) 
and the initial conditions: at t = 0, 


(4.1) 



if x > 0 

if x < 0 


wiiere U and V are two different given constants. For t. > 0, the exact weak solution 
to this problem is 


(4.2) 


f U, if x — at > 0 
\ V, if x — at < 0 


Obviously, the space-time variation of the above solution is solely driven by the 
initial-data discontinuity that occurs at x = 0. 

To construct a corresponding discretized initial- value problem, consider an un- 
bounded and uniform space-time mesh formed by the mesh points (j, n) £ ft with 
n > 0 (see Fig. 6). The numerical analogue of the initial conditions Eq. (4.1) can 
be expressed as (i) 


(4.3) 

and (ii) 

(4.4) 


7( o_ J U, if j = 1/2, 3/2,... 

3 \ V, if j = -1/2, -3/2,... 

(«+)? = 0, for j = ± 1/2, ±3/2,.... 


Let ft+ be the subset of ft with n > 0. Then, for any (j, n) £ ft+, u ” and ( u +)” 
can be determined in terms of the above given initial data through the marching of 
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a given a' scheme. Hereafter the solution thus obtained and the associated initial- 
value problem will be referred to as Solution I and Problem I, respectively. As in 
the case of its analytical counterpart, the space- time variation of Solution I is also 
solely driven by the initial-data discontinuity that occurs across the mesh interval 
centered at x = 0. 

The mesh depicted in Fig. 6 is unbounded. In reality one can only use a 
bounded mesh. To proceed, first we shall define several finite subsets of f7 + . Let 
(i) j 0 > I and n 0 > 0 be given whole integers, and (ii) j b c = f j 0 ± 1/2. Then the 
mesh-point set A+(jo,no) and its subsets A' + (jo,no) and <9A+(jo,no) are defined as 
follows: (i) (j,n) E A+(j 0 ,n 0 ) if and only if (j,n) E Q+, |j| < j b and n < n 0 4- 1/2 
(i.e., A + (jo,no) is the set of the mesh points (j,n) depicted in Fig. 7 excluding 
those with n = 0); (ii) (j,n) E A+(jo,no) if and only if (j, n) E A+(jo,no) and 
| j\ < jb , and (iii) (j,n) E 3A + (j 0 ,n 0 ) if and only if (j,n) E A+(j 0 ,n 0 ) and \j\ = j b . 
Obviously A + (jo,n 0 ) = A^(j 0 ,n 0 ) U 9A + (j 0 ,ri 0 ). 

Next we cosider the initial data defined by (i) 


Note that u j and (u+)j defined in Eqs. (4.5) and (4.6) are identical to those defined 
in Eqs. (4.3) and (4.4), respectively, in their common domain defined by |j| < j b . 
The question arises that, given only the truncated initial data Eqs. (4.5) and (4.6), 
is it still possible to obtain a solution in A' + (jo,rio) that is reasonably close to 
Solution I? 

To answer the above question, first consider a discrete initial/boundary- value 
problem in which the initial data are specified using Eqs. (4.5) and (4.6), and the 
boundary data a” and ('u+)™ at all (j, n) E 5A + (jo,no) are also given. Because, for 
any (j, n) E 17, (j, ra) and (j ± 1/2, n — 1/2) form the stencil of the a ' scheme, one 
concludes that, for all (j,n) E A+(j 0 ,no), u™ and (aj)t' can be uniquely determined 
in terms of the above given initial and boundary data by using the same a! scheme 
that was used to generate Solution I. Hereafter, the new solution thus obtained 
and the associated initial/boundary-value problem will be referred to as Solution II 
and Problem II, respectively. Obviously, Solution II is dependent on the boundary 
data specified on 3A+(j 0 ,no)- By definition, the degree of non-reflectiveness of 
the boundary conditions of Problem II will be measured by how closely Solutions 
I and II are matched in A+(jo,no). In other words, the closer the match the less 
reflective these boundary conditions are. In particular, they are said to be perfectly 
nonreflecting if the match is perfect. 

Next note that, again because of the shape of the stencil of an a ' scheme de- 
scribed earlier, for Solution I, u 7 - and ( u +)? at all (j, n) E A+(jo,n 0 ) are completed 
determined by the initial data given at the mesh points (j, 0), j = ±1/2, ±3/2, . . . , 
±(j b + n 0 ) (see Figs. 6 and 7). As a result, one can see easily that Solutions I 
and II must be identical in A' + (jo,no) if the values of Solution I at 3A + (jo,no) are 
used as the boundary data for Problem II. In other words, the boundary conditions 
formed using these boundary data are perfectly nonreflecting for Problem II. Un- 
fortunately, obtaining these boundary data requires solving a problem with a larger 
spatial domain, whose size increases as n 0 increases. 


(4.5) 



and (ii) 
(4.6) 


(05 = 0, for j = ±1/2, ±3/2, . . . , ±jb 
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The key propositions that, emerge from the above discussions are: 

(a) Any set of NRBCs for Problem II must meet the requirement that the re- 
sulting Solution II more or less matches Solution I in A' + (jo,n 0 ). 

(b) In constructing any set of NRBCs for Problem II, only the information 
extracted from the initial data Eqs. (4.5) and (4.G) can be used. 

(c) The values of Solution I in A' + (jo, no) generally are dependent on some initial 
data beyond those given in Eqs. (4.5) and (4.6). 

It follows from these propositions that construction of NRBCs for Problem II is 
possible only if the contribution to Solution I in A+(jo,no) by the initial data Uj 
and (u+)j with \j\ > jb somehow can be extrapolated based on the truncated 
initial data Eq. (4.5) and (4.6). Fortunately, the existence of such an extrapolation 
is supported by the following facts: 

(a) Except in the neighborhood \j\ < 1/2, the initial data Eqs. (4.5) and (4.6) do 
not vary in each of the remaining regions 1/2 < j < jb and —jb < j < — 1/2. 

(b) Not only does the initial data Eqs. (4.3) and (4.4) coincide with those of 
Eqs. (4.5) and (4.6) in their common region \j\ < jb , in the region j > Jb 
(j < — jb ), the former initial data also have the same constant values of those 
in the region 1/2 < j < j b (-jb < j < -1/2/. 

In the following, it will be shown that practical, albeit not perfect, nonreflecting 
boundary conditions for Problem II can be constructed using an assumption about 
the flux distribution near its spatial boundaries. The exact meaning of this assump- 
tion will become completely clear only after its role in constructing NRBCs is shown 
explicitly. However the gist of this assumption will be described immediately. 

As an example, consider the mesh line j = jb, i.e., the right boundary mesh line 
of Problem II depicted in Figs. 6 and 7. Given Solution I, one can evaluate, from 
t = 0 to t. = rioAt. the distribution of the space-time flux of h* passing through (say 
the positive direction is from left to right) the above mesh line and its neighboring 
vertical mesh lines. Because the space-time variation of Solution I is solely driven 
by an initial-data discontinuity that is located to the left of the mesh line j = jb, 
no propagating “solution disturbances” will ever reach this mesh line from its right 
side. As a result, one may expect that, for the time interval (0,rtoAt), the flux 
distribution of Solution I at the mesh line j = jb is not “dependent” on those at 
the neighboring vertical mesh lines with j > jb, i.e., it is “ dependent ” only on those 
at the neighboring vertical mesh lines with j < jb. In fact one may further assume 
that, for the time interval (0, no At), the flux distribution of Solution I at the mesh 
line j = jb is a smooth extrapolation of those at the neighboring mesh lines with 
j < jb. Because a set of NRBCs for Problem II is introduced such that the resulting 
Solution II more or less matches Solution I in A' + (jo, no), it follows that Solution II 
must possess the same flux extrapolation relation referred to above if the boundary 
conditions imposed at the mesh line j = jb are indeed nonreflecting. As will be 
shown shortly, the current NRBCs are derived assuming this flux extrapolation 
relation. Note that, according to the current line of thinking, a set of NRBCs is 
one that allows for the flux to smoothly “stream” out of the spatial boundary from 
the interior. 

To elaborate further, consider a ID problem (-boo > x > — oo) that is governed 
by Eq. (2.1). For such a problem, solution disturbances propagate to the right 
(left) with a constant speed if a > 0 (a < 0). As a result, no solution disturbances 
will ever reach the mesh line j = jb from its right side if the initial sources of 
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disturbance are located to the left of the mesh line. In contrast, for a more complex 
problem such as that involving a ID subsonic Euler flow, (i) disturbances emitted 
from any source, relative to the stationary computational frame , can propagate in 
both the right and left directions; and (ii) any disturbance itself is a new source 
of disturbance. As a result, for such a flow, solution disturbances that originate 
from a source that is located to the left of the mesh line j = ji> can propagate into 
the region to the right of the mesh line and from there propagate back to the same 
mesh line. In the current construction of NRBCs, the effect of this “back scattering” 
phenomenon on the flux distribution along the mesh line j = jt, is assumed to be 
negligible. Generally speaking, the last assumption is valid if the mesh point (0, jb) 
is located at a large distance from the nearest initial source of disturbance (such 
as in the current case in which by assumption > 1, and the initial source of 
disturbance is located around j = 0). 

This subsection is concluded with the following remarks: 

(a) The above discussions make it clear that, strictly speaking, NRBCs for a 
discretized initial/boundary- value problem are not well-defined unless its 
solution can be compared against that of another problem with a space- 
time domain containing that of the original problem. 

(b) Conceptually, the reflecting boundary conditions that are imposed on a solid 
wall are exactly opposite to the NRBCs discussed above. Obviously, the 
mass flux cannot smoothly “stream” into a solid wall. 

(c) In the following subsections, the concept of GCEs introduced in Sec. 3 will 
be used in the construction of several sets of NR BCs which are applicable to 
any a' scheme. Before proceeding to construct these conditions, the reader 
again is reminded of the convention that any part of the boundary of a GCE 
on which h* is not specified at each individual point is denoted by a dashed 
line (see Fig. 4). 

4.2. The First and Second Sets of NRBCs. The boundary mesh lines 
with j = ±jb that appear in Figs. 6 and 7 along with some of their immediately 
adjacent interior mesh lines are depicted in Fig. 8. In this figure, the points of 
intersection between these vertical mesh lines and any horizontal nth time-level 
line, respectively, are denoted by At, Bi,Ct,Di,A! t , B'^C[ and D \ , where £ — 2 n. 
Furthermore, for each l = 0, 1, 2, . . . , points Pi, Qi, P[ and Q ' t , respectively, are on 
the line segments BfC( , C(D( , B[C[ and C'^D^ with 

(4.7) 

\m\ = ^, |cw7l = (i-A)^, \P' l c' l \ = x'^, |qe'| = (i-A')^ 

Here A and A r are adjustable parameters with 0 < A, A' < 1. Obviously \PiQe\ = 
\PfQW = ax/ 2. Furthermore, because \BiCi\ = \CiDi\ = ax/2, one concludes 
that, for i — 0,1,2,..., points Pi and Qi , respectively, coincide with (i) points 
Bf and Ci if A = 1 , and (ii) points Ce and Di if A = 0. As a result, as long as 
0 < A < 1, for each i — 0,2,4,..., the interior of PiQt Ues entirely within the 
solution element associated with point Ci. Obviously the above conclusions also 
hold if points Pi , Qi , Bi, Ci and D ( , and the parameters A are replaced by points 
P' ( , Q' t , B\ , C[ and D\, and the parameter A', respectively. 

According to the discussions given in Sec. 3, for any £ = 0, 1, 2, . . . , the rectan- 
gles AiAi+\ Bi+\ Bi , PiPt+\ Ct+\ Ci, CiCi+\ Qi+\ Qi , B' f , P[P{+\ C /+ 1 C[ 
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and can bo considered as GCEs. As a result, the union of any com- 

binations of these GCEs is also a GCE. As an example, for any £ = 0,2,4,..., 
the rectangle AiAt+ 2 Bf+ 2 B( is a GCE. To proceed, for any £ = 0,2,4,..., let (i) 

F(A e A w ) denote the flux of h* entering the vertical line segment AiAi+2 from 
the left , and (ii) F(AfBf) denote the flux of h * entering the horizontal line segment 
AfBf from below. Hereafter, similar notations will be used to denote the fluxes 
passing through other vertical and horizontal line segments. The simplest NRBCs 

at the boundary mesh line j = jt> will be constructed by assuming the simplest flux 
extrapolation condition, i.e., 

(4.8) £ = 0, 2, 4, . . . 

Because the rectangle AiAt+ 2 Bt+ 2 Bi is a GCE, we have 

(4.9) F(AeA w ) ~ F(B,B w ) + F{A M B t + 2 ) - F(AjB~ f ) = 0, £ = 0, 2, 4, ... . 
Combining Eqs. (4.8) and (4.9), one concludes that 

(4.10) F(A f+2 B (+2 ) = F(A&i), £ = 0, 2, 4, . . . 

As a result, 

(4.11) F(AtB f ) = F(OTo), £ = 2, 4, 6, . . . 

With the aid of Eqs. (2.4), (2.5) and (2.9), and recalling that the length of each 
AfBf is ax/ 2, it can be shown that Eq. (4.11) is equivalent to 

(4.12) u Ae ~(u+) Ae =UA 0 ~(u+)ao, £ = 2,4,6... 

where ua, and (u+)A f , £ = 0, 2, 4, . . . , denote the mesh values of u and u+ at point 
Af, respectively. Note that: (i) hereafter similar notations will be used to denote 
the mesh values of u and u+ at other mesh points; and (ii) according to Eq. (2.9), 
( u t ) a ( and {u+)a 0 are quantities of first-order in ax. Also note that, by using a 
similar approach, it can be shown that the counterpart to Eq. (4.12) for the left 
boundary is 

(4.13) u A ’ t + {u+) A > t = u a > q + (u$) A ' 0 , £ = 2,4,6,... 

Let the initial data ila 0 and (u+)a 0 given. Then, for each £, there are only 
two unknowns, i.e., the boundary data ua ( and (uj)^, in Eq. (4.12). As a result, 
after imposing Eq. (4.12), there is still one degree of freedom left in choosing the 
two unknowns. To obtain the simplest boundary conditions, the degree of freedom 
is removed by further assuming that the zero-order term and the first-order term 
on the left side of the equality sign in Eq. (4.12) are equal to those on the right 
side, respectively. Thus one has 

(4.14) u At =u A 0 arid (u+) Al = (u+) Ao £ = 2,4,6,... 

Similarly, the counterpart to Eq. (4.14) for the left boundary is 

(4.15) u A ’ t =u A ' o and («+) A j = (u+) A > 0 £ = 2,4,6... 

Note that, even though it is counter-intuitive, it will be shown in Sec. 5 that the 
Euler counterparts of the steady-state boundary conditions Eqs. (4.14) and (4.15) 
turn out to be a set of reasonably good NRBCs. Hereafter Eqs. (4.14) and (4.15) 
will be referred to as the first set of NRBCs. 
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A set of more advanced NRBCs will be constructed assuming 


(4.10) F(AtAt+i) - F{D t D w ) = F{P e P W ) - F(Q ( Q e + 2 ), t = 0, 2, 4, . . . 

Note that, because (i) the line segments A[Ai+ 2 , D( £^+ 2 , P(Pe+ 2 ? and QfQf-f 2 have 
the same length (i.e., At), and (ii) the horizontal distance (i.e., ax/2) that separates 
the line segment AtAt+ 2 and B(B (+2 is identical to that which separates the line 
segments P 1 P 1+2 and QtQt + 2 , Eq. (4.10) represents a linear flux extrapolation 
relation. Also, because the rectangles AtAt+zBe+zBf and PfPt+iQf+iQi are GCEs, 
we have Eq. (4.9) and 

(4.17) F(P^T 2 ) - F(QjQ^ 2 ) + F(P e +2Qt+2) - F(7W() = 0 

for t — 0,2,4, — With the aid of Eqs. (4.9) and (4.17), Eq. (4.10) implies that 

(4.18) F(A t+2 B w ) - F{P ( + 2 Q(+2) = F(AiBi) - F{T\Qc), l = 0, 2, 4, . . . 

As a result, 

(4.19) F(MB ( ) - F(TW<) = F(OTb) - F(P^), t = 2, 4, 6, . . . 

With the aid of Eqs. (2.4), (2.5), (2.9) and (4.7), it can be shown that Eq. (4.19) is 
equivalent to 

(4.20) (u A , - u Ct ) - [(ut)A' + (2A - l)(u+) c< ] 

= (UA 0 ~ UC 0 ) ~ [(«J)^o + ( 2A “ W^Co] 

where ( = 2,4,0, Similarly, the counterpart to Eq. (4.20) for the left boundary 

is 

(4.21) (u A , - u c .) + [(«+)*; + (2A' - l)(«+)cj] 

= ( U A‘ 0 - + [(uJJaj, + (2A' - l)(w^)cj 

where t = 2,4,0, — To proceed, note that Eqs. (4.5) and (4.0) coupled with the 
assumption > 1 imply that ua 0 = uc 0 = U , ua> 0 = uc' 0 = V and (uj)^ 0 = 
(^)co = ( u t)A' 0 = (ti+)cj =0- As a result, the expressions on the right sides of 
Eqs. (4.20) and (4.21) vanish. 

Let the parameter A be given. Then because the Ct is an interior mesh point 
at the (^/2)th time level (recall that n = 1/ 2), uc t and (u+)c e can be evaluated 
using the known marching variables at the ((f. - l)/2)th time level. Thus one 
concludes that, for each t, there are only two unknowns, i.e., the boundary data 
UA f and (;Uj )a,. in Eq. (4.20). As a result, after imposing Eq. (4.20), there is still 
one degree of freedom left in choosing the two unknowns. To obtain the simplest 
NRBCs possible, the degree of freedom is removed by further assuming that the 
sum of the zero-order terms and that of the first-order terms on the left side of 
Eq. (4.20), respectively, are equal to zero. Thus one has 

(4.22) u Al = uc t and (ut) A , = (1 - 2A)(u+)c, (= 2,4,6,... 

Similarly, the counterpart to Eq. (4.22) for the left boundary is 

(4.23) u A ' t =u C ' t and (u+) A - ( = (1 - 2\')(u+)c t ^ = 2,4,6,... 

Note that, because 1 > A, A' > 0, 

(4.24) 1 > (1 - 2A), (1 - 2A') > -1 

Hereafter Eqs. (4.22) and (4.23) will be referred to as the second set of NRBCs. 
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4.3. The Third Set of NRBCs. The points Af , Bf, Cf , A ' ( , and C^, 

£ = 0, 1, 2 , that are depicted in Fig. 9 are identical to those depicted in Fig. 8. 

Furthermore, for each £ = 0, 1, 2, ... , points Rf , Sf , and S' ( , respectively, are on 
the line segments AfBf, BfCf , A\.B\ and B' f C' ( with 

(4.25) 

At at At At 

\R ( B(\ = X—, |B3£| = (1-A) t . |^| = A' t , |B^| = (1 -A') t 

Here A and A' are adjustable parameters with 0 < A, A' < 1. Obviously \RfSf\ = 
I^S^I = ax/2. Furthermore, because \A f Bf\ = = at/ 2, one concludes that, 

for £ = 0, 1, 2, . . . , points Rf and Sf , respectively, coincide with (i) points Af and 
B( if A = 1, and (ii) points Bf and Cf if A = 0. As a result, as long as 0 < A < 1, 
for each £ — 1, 3, 5, ... , the interior of RfSf lies entirely within the solution element 
associated with point Bf. Obviously the above conclusions also hold if points Rf , 
Sf , Af , B f and Cf , and the parameters A are replaced by points R' f , Sf, A! t , B\ and 
C’ f , and the parameter A', respectively. 

The third set of NRBCs will be constructed assuming 

(4.26) F{A/Aj^) - F(B/T3^ 2 ) = F(R^~JT^) - F{SiI^S/^), £-2,4,6,... 

Note that, because (i) the line segments AfAf+ 2 * BiBe+ 2 , Rf-iRi+i , and Sf-iSf+i 
have the same length (i.e., At), and (ii) the horizontal distance (i.e., at/ 2) that sep- 
arates the line segment AfAf + 2 and BfBf + 2 is identical to that which separates the 
line segments Rf^iRf+i and Sf-iSf+\, Eq. (4.26) represents a linear flux extrapo- 
lation relation. Also because the rectangles AfAf^Bi^Bf and Rf-iRf+iSf+iSf-i 
are GCEs, we have Eq. (4.9) and 

(4.27) F(Rf-\Rf+\) - F(Sf-\Sf+\) + F(R W S W ) - F(R^I 7) = 0 

for £ = 2,4,6, With the aid of Eqs. (4.9) and (4.27), Eq. (4.26) implies that 

(4.28) 

F(Af +2 Bf+ 2 ) — F(Re+iSi+i) = F(AfBf) - F(Ri_iSf- 1), £ — 2,4,6, . . . 

As a result, 

(4.29) F(A^Bl) - F(Rf-iSf-i) = F(MB~ 2 ) - F(R ^), £ = 4, 6, 8, . . . 

With the aid of Eqs. (2.4), (2.5), (2.9), and (4.25), it can be shown that Eq. (4.29) 
is equivalent to 

(4.30) (u At - «/?,_,) - + (2A - l)(u+) B ,_,] = (ua 2 -u B i ) 

- [{u+)a 2 + (2A - 1 )(«+)* ] , e = 4, 6, 8, . . . 

Similarity, the counterpart to Eq. (4.30) for the left, boundary is 

(4.31) (ua 1 , - «();_,) + [( w J)4i + (2A' - l)(u+ )c;_J = (u^ — «b;) 

+ [(<4Ui + (2A'-l)(«+)*j], / = 4,6,8,... 

To proceed, note that: (i) Eqs. (4.3) and (4.4) coupled with the assumption jb 1 
imply that, for Solution I, 

(4.32) u°j = U and (u+)° =0, j = jbji, ± 1 
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and (ii) each a ' scheme has the property that, for any ( j, n) 6 ft, u “ = and 

(u+)” = 0 if y^Zi /2 = w j+i /2 an d ( w x )j±i 1 /2 = 0* As a result, for Solution I, 

(4.33) »)Zm = U “ d =° 

and 

(4.34) < = U and (w+)j 6 =0 

Because, in Figs. 8 and 9, the mesh points ( jt, — 1/2, 1/2) and (j&,l) are denoted 
by D\ and respectively, Eqs. (4.33) and (4.34) imply that, for Solution I, 

(4.35) im, = tig, = U and (u+)a 2 = (u+) B , = 0 
Similarly, it can be shown that, for Solution I, 

(4.36) u A ’ 2 = u B \=V and (u+) A > 2 = (u+) B ; = 0 

It will be assumed that Eqs. (4.35) and (4.36) are also valid for Problem II. As a 
result, the expressions on the right sides of Eqs. (4.30) and (4.31) vanish. 

Let the parameter A and the marching variables ub £-1 and (u+)b ( _ 1 (which 
are associated with a lower time level) in Eq. (4.30) be given. Then, for each 
f. = 4,6,8,..., there are only two unknowns, i.e., the boundary data ua, and 
( u t ) a ( in Eq. (4.30). As a result, after imposing Eq. (4.30), there is still one degree 
of freedom left in choosing the two unknowns. To obtain the simplest NRBCs, the 
degree of freedom is removed by further assuming that the sum of the zero-order 
terms and that of the first-order terms on the left side of Eq. (4.30), respect uively, 
are equal to zero. Thus one has 

(4.37) u Al =u Bl _ 1 and (w+) 4< = (1 - 2A)(uJ) c< _ 1 1 = 2,4,6,... 

Similarly, the counterpart of Eq. (4.37) for the left boundary is 

(4.38) u A ' t =u B ' t _ , and («+)*»= (1 - 2 A')(u+)b'_, ^ = 2,4,6,... 

Note that: (i) obviously Eq. (4.24) must also be observed here; and (ii) hereafter 
Eqs. (4.37) and (4.38) shall be designated as the third set of NRBCs. 

4.4. The Fourth Set of NRBCs. Each of the three sets of NRBCs intro- 
duced above represents only the simplest among many possible combinations of 
boundary conditions that are consistent w r ith the associated flux extrapolation con- 
ditions. As an example, instead of using the simplest approach described in Sec. 4.2, 
nonuniqueness involving Eqs. (4.20) and (4.21) can also be removed by assuming 
that, for I = 2,4,6 , 

(4.39) UA t — uc t -f a x(u x )c t and ua' ( = Uc t — Ax(u x )c l— 2, 4, 6, . . . 

i.e., ua, and UA' t are related to uc t and uc' e , respectively, by the first-order Tay- 
lor’s expansion. Combining Eqs. (4.20), (4.21), and (4.39), and the fact that the 
expressions on the right sides of Eqs. (4.20) and (4.21) vanish, it can be shown that 

(4.40) 

(ut) At = (5 - 2A)(u+)c, and (u+) A j = (5 - 2A')(u+)cj, / = 2,4,6,... 

Eqs. (4.39) and (4.40) shall be designated as the fourth set of NRBCs. Note that, 
because 1 > A, A' > 0, one has 

(4.41) 5 > (5 — 2A), (5 — 2A # ) >3 


NAS A/TM— 2003-2 1 2495/REV 1 


15 


4.5. Remarks and Extensions. This section is concluded with the following 
comments: 

(a) To simplify discussion, it is assumed in Eqs. (4.3)-(4.6) that the initial source 
of disturbance is a single discontinuity which is at a reasonable distance away 
from the mesh points (0, ±j&). However, according to the logic of the current 
development, the NRBCs derived here are still valid even if the discontinuity 
is replaced by a source region of finite extent as long as this source region is 
kept at a reasonable distance away from any nonreflecting boundary. 

(b) With the aid of Eqs. (2.11) (2.15) and the related remarks, the develop- 
ment presented in this section can be extended to the ID Euler case in a 
straightforward manner. Specifically, in the Euler development, (i) Eq. (2.1) 
is replaced by Eq. (2.11); (ii) Eqs. (4.1) and (4.3)-(4.6), respectively, are re- 
placed by their Euler versions that, respectively, are identical to these equa- 
tions except that u, u°, (u+)j, U, and V respectively, are replaced by u m , 
( u m )j , (u+Jj, U m , and V m , where m = 1,2,3, and U m and V m are two 
different constants for each m; (iii) the exact solution given in Eq. (4.2) is 
replaced by its Euler version, i.e., the Riemann solution [41, p. 181]; and 
(iv) the condition 0 < A, A' < 1 is again assumed. Using the new definitions 
given in the above item (ii), one can proceed to define the Euler versions of 
Problems I and II, and Solutions I and II. Obviously, (i) by replacing the 
symbols u and with ( u m ) and (u+ x ), respectively, the NRBCs and (u+ x ), 
respectively, the NRBCs Eqs. (4.12) (4.15), (4.20) (4.23), (4.30), (4.31), and 
(4.37) (4.40) will become their Euler versions, respectively. Note that here- 
after the Euler version of an equation such as Eq. (4.3) will be denoted by 
Eq. (E-4.3). 

(c) Consider the ID Euler case. Because m = 1,2,3, each of Eqs. (E-4.12), 
(E-4.13), (E-4.20), (E-4.21), (E-4.30), and (E-4.31) represents a set of three 
conditions involving six unknowns. Thus, after imposing any one of these 
equations at a boundary mesh point (j,n), there are still three degrees left in 
specifying (u m )y and (u+ x )™. The several approaches used here to eliminate 
these degrees of freedom represent only a small subset of all possible options. 
It is conceivable that, by using other options, one may be able to construct 
a set of NRBCs that also meet some extra given conditions at the boundary. 

5. Numerical Results 

In this section, the Euler versions of the four sets of NRBCs developed in Sec. 4 
will be validated numerically. Because only the Euler versions are considered here, 
in this section a phrase such as “the Euler version of Problem I” will be abbreviated 
simply as 1 Problem I”. 

The ID test problem to be specified shortly is an extended Sod’s shock tube 
problem. It is the original Sod’s problem [42] with the additional complication of 
imposing NRBCs at the two open ends of the shock tube. Note that the flow under 
consideration contains a shock wave and a contact discontinuity, and, relative to 
the computational frame , is subsonic throughout (the values of Mach number range 
from 0 to 0.93). It is well known that implementing NRBCs for a subsonic flow is 
much more difficult than doing the same for a supersonic flow. This difficulty is 
further exacerbated by the fact that most established NRBCs are derived assuming 
a continuous flow — which is not valid for the current case. In this section, each of 
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the four sets of NR DCs developed in Section 4 will be validated by comparing the 
resulting numerical solutions against the exact analytical solutions for a period of 
time that ends only after the shock wave and contact discontinuity have exited the 
computational domain. 

To proceed, let the meshes used by Problems I and II, respectively, be those 
depicted in Figs. 6 and 7 with a t = 0.004, ax = 0.01 and jo = 50. Thus, for 
Problem II, jt, = 50| and —0.505 < x < 0.505. In addition, it is assumed that: 
(i) the specific heat ratio 7 = 1.4; (ii) both Solutions I and II are obtained using 
the ID Euler a-a scheme with a = 1; and (iii) the initial conditions at t = 0 for 
Problem I (Problem II) are Eqs. (E-4.3) and (E-4.4) (Eqs. (E-4.5) and (E-4.G)) with 
(UuU 2 ,U 3 ) = (0.125, 0, 0.25) and (V U V 2 ,V 3 ) = (1, 0, 2.5). For both Solutions I 
and II, the above assumptions imply a Courant number CFL ~ 0.88 where CFL is 
defined to be the maximal numerical value of of (|u| 4- c)At/Ax with v and c being 
the local velocity and sonic speed, respectively. Note that stability of the ID Euler 
a-a scheme generally requires CFL < 1 and a > 0. 

Moreover, by using (i) Fig. 6; (ii) Eqs. (E-4.3) and (E-4.4); and (iii) a property 
of the ID Euler a-a scheme, i.e., for any ( j,n ) E D, (u rn )^ = (ttm)”^ 2 and 

WJj = 0 if {u m )”Zl/2 = 311(1 («m*)"±l/2 = °> jt Can be sh ° WI1 that ’ 

for Solution I, 


(5.1) 

and 

(Um)? = U m 

and 

(«m*)" = 0 

if j > n 

(5.2) 

( U m) r j = 

and 

( u mx)j = 0 

if j < -n 


In other words, the impact of the initial discontinuity is felt only in the mesh region 
\j\ < n (note: \j\ ^ n if (j,n) E D). Because jt, = j 0 + 1/2, Eqs. (5.1) and (5.2) 
imply that 

(5.3) {umYj b = Um and (u+ x )" b = 0, n = 0, 1,2, . . . , j 0 
and 

(5.4) {u m )-j h = Vm and (u+ x )" jb =0, n = 0, 1,2 ,.. . ,j 0 

respectively. Furthermore, because (i) n = ^/2, (ii) (u m )A 0 = ( u m)j k = U m 
and (u+ x ) Ao = (“!«)“, = 0, and (iii) (u m )A' 0 = (u m )- jb = V m and (u+ x ) A ' 0 = 
( u mi)-jfc = 0, a comparison of Eqs. (E-4.14) and (E-4.15) with Eqs. (5.3) and (5.4) 
reveals that, for n = 0, 1, 2 , . . . , jo, the boundary values at the mesh lines j = ±jb 
specified using the first set of NRBCs are identical to those specified using the so- 
lution values of Solution I. According to a discussion given in Sec. 4.1, this implies 
that Solution II coincides with Solution I at all mesh points (j, n) with |j| < jt, and 
0 < n < jo if the boundary values of Problem II are specified using the first set of 
NRBCs. Note that, given any pair of the values of A and A', one can also show that: 
(i) Solution II coincides with Solution I at all mesh points (j, n) with |j| < jt, and 
0 < n < jo if the third set of NRBCs is used; and (ii) excluding the two boundary 
mesh points with j = ±jb and n = jo, Solution II coincides with Solution I at all 
mesh points (j, n) with |j| < jt, arid 0 < n < jo if the second or the fourth set of 
NRBCs is used. 

Because jo = 50 and At = 0.004, t = nAt = 0.2 when n = jo- As a result, 
the above discussions imply that, within their common spatial domain, Solution 
II should be identical to Solution I in the time interval 0 < t < 0.2 if the first or 


NASA/TM— 2003-2 1 2495/REV 1 


17 


the third set of NRBCs is used. Also, except for the two boundary mesh points 
mentioned above, the above conclusion also applies if the second or the fourth set 
of NRBCs is used. Note that, in reality, the differences between Solutions I and II 
are completely negligible at the two exceptional boundary mesh points. 

At t = 0.2, Solution II generated using the first set of NRBCs is shown in 
Fig. 10(a). In this figure, the numerical values (marked by triangles) of the density, 
the velocity and the pressure are compared with the exact solutions (marked by 
solid lines). As expected from the discussions given in the last paragraph, at t — 
0.2, Solution II generated using any one of the other three sets of NRBCs is also 
represented by the same results shown in Fig. 10(a). Note that, at t. = 0.2, the waves 
and shocks generated in the interior have not yet reached the boundaries. Also it 
is seen that the agreement between the numerical results and the exact solutions is 
excellent. In particular, the shock discontinuity is resolved almost within one mesh 
interval, and the contact discontinuity is resolved in four mesh intervals. Also, there 
are only slight numerical overshoots and/or oscillations near these discontinuities. 

At t — 0.4, Solution II generated using the first set of NRBCs is shown in 
Fig. 10(b). It is seen that, by this time, the shock wave has passed cleanly through 
the right boundary. There is good agreement between the numerical solutions and 
the exact solutions everywhere in the interior except for a slight disagreement in 
the vicinity of the right boundary. Note that the right boundary values , which do 
not vary with time and are no longer identical to the values of Solution I there 
at t = 0.4, are discontinuous with respect to the neighboring interior values. The 
numerical results at t = 0.6 are shown in Fig. 10(c). As seen from the density 
profile, by this time, the contact discontinuity has also passed through the right 
boundary. Agreement between the numerical solutions and the analytical solutions 
continues to be good in the interior. However, both left and right boundary values 
are now discontinuous with respect to the neighboring interior values. 

At t = 0.2, 0.4, 0.6, solution II generated using the third set of NRBCs with A = 
A' = 0 are shown in Figs. 10(a), 11(a) and 11(b), respectively. The agreement with 
the exact solution is excellent everywhere. In particular, unlike the case associated 
with the cruder first set of NRBCs, agreement between the numerical and exact 
solutions at t = 0.4 and t = 0.6 are quite good for the current case even at the right 
boundary. Note that the use of other values of A and A' in the range 0 < A, A' < 1 
yields almost identical numerical results. 

To give the reader a concrete idea about the simplicity of the CE/SE method in 
general and the current NRBCs in particular, the self-contained Fortran program 
which is used to generate the results depicted in Figs. 10(a), 11(a) and 11(b) is 
listed in Appendix. In this program, the input and output include: (i) nx = 2 jb 
= the number of total mesh intervals at the 0th time level, (ii) it = the number of 
total marching steps (each marching step advances a time period of At/ 2), (iii) dt 
= At , (iv) dx = ax, (v) (ul,u2,u3) = (U \ , t/ 2 , t/ 3 ), (vi) (vl,v2,v3) = (V i,V^V 3 ), 
(vii) cr = 1 — 2A, (viii) cl = 1 — 2A', (ix) ia = a, (x) rho = mass density, (xi) 
v = fluid velocity, (xii) p = static pressure, and (xiii) M = Mach number. Note 
that (1 - 2A) and (1 - 2A') are the factors that appear in Eqs. (4.37) and (4.38), 
respectively. Thus, according to Eq. (4.24), 1 > cr, cl > —1. 

Solution II also has been generated using the second and the fourth sets of 
NRBCs with different values of A and A'. At t = 0.2, 0.4, 0.6, the numerical results 
thus obtained are virtually indistingushible (at least visually) from those shown in 
Fig. 10(a), 11(a) and 11(b) for all A and A' in the range 0 < A, A' < 1. 
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Note that the second, the third and the fourth sets of NRBCs were derived 
assuming 0 < A, A' < 1. As a result, breakdown of these NRBCs is expected if the 
values of A and A' used fall outside of the range 0 < A, A' < 1. As an example, 
consider the fourth set of NRBCs with A = A' = 2. For this special case, Eqs. (E- 
4.39) and (E-4.40) can be rearranged as 

(5.5) (UrnjAt = ( u m)ce 4" ^ x { u mx)Ce ( u mx) At — ( u mx)Ci 

and 

(5.6) (Um)^ = ( u m)cg ~ Ax(n m x)cj ( u rnx)A' g = ( u mx)c , ( 

where m = 1,2,3 and t = 2,4,6, Because (i) (u m )A e and (u m )A' f are related 

to (u m )c e and (u m )c' e , respectively, by the first-order Taylor’s expansion; and (ii) 
{ u mx)^e ( u mx)A’ ( simply assume the values of the corresponding variables at 
the neighboring mesh points Cf and C[ y respectively, these boundary conditions 
certainly look “reasonable”. As explained earlier, at t = 0.2, Solution II gener- 
ated using Eqs. (5.5) and (5.6) is again represented by those shown in Fig. 10(a). 
However, at t = 0.4 (see Fig. 12), Solution II generated using these “reasonable” 
boundary conditions has become highly reflecting at the right boundary. 

6. Conclusions and Discussions 

In this paper, we first review the CE/SE schemes for the ID advection equation 
Eq. (2.1) and those for the ID Euler equations Eq. (2.11). The concept of a general- 
ized conservation element was then developed, paving the way for the development 
of flux-based NRBCs. 

Four sets of simple NRBCs for the advection equation were developed on the 
basis of extrapolation of fluxes at the boundaries. No characteristics-based tech- 
niques were used in the development. The Euler versions of these NRBCs w r ere 
developed in a similar manner and tested with the ID Euler a-a scheme. The test 
problem used is an extension of the original Sod’s shock- tube problem in which 
extra NRBCs are imposed at the two open ends of the shock tube. 

The first set of Euler NRBCs is the simplest and represents a set of steady- 
state boundary conditions. It was truly remarkable that this set of NRBCs yielded 
numerical solutions that were in good agreement with the exact w r eak solution in 
the interior of the computational domain , even after the prescribed steady-state 
boundary values had deviated completely from those associated with the exact 
solution. 

It was also shown that the more advanced second, third and fourth sets of Euler 
NRBCs yielded numerical solutions that were in even better agreement with the 
exact weak solution — if the values of the adjustable parameters A and A' used in 
these NRBCs fell within the allowable range 0 < A, A' < 1. However, as shown by 
a numerical example, strong spurious reflections from a boundary may occur if the 
values of A and A' chosen fall outside the allowable range. 

Finally note that it has been shown numerically [19,21 28] that the multidi- 
mensional extensions of the current Euler NRBCs are also equally simple, effective 
and robust. In particular, the reader is referred to [19] for a rigorous discussion 
of the investigation of several classes of flow problems using a 2D CE/SE Euler 
solver [15] in conjunction with a 2D extension of the present ID Euler NRBCs. 
The classes of problems considered include (i) acoustic pulse, entropy wave, and 
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vorticity wave propagation, (ii) free shear-layer instability, and (iii) multiple inter- 
actions of a strong vortex and strong oblique shocks. The major conclusions of this 
investigation are: 

(a) The (second-order) 2D CE/SE solver is efficient and yields high resolution, 
low dispersion results similar to those of higher-order conventional schemes 

[4-6]. 

(b) The novel 2D CE/SE NRBCs used are easy to implement and remove the 
need for (or mimimize the size of) a buffer zone. 

(c) The 2D CE/SE solver is capable of handling both continuous and discon- 
tinuous flows and, thus, provides a unique numerical tool for flow situations 
where sound waves and shocks and their interactions are important, such as 
the jet screech noise problem. 

Appendix: A Solver for An Extended Sod’s Shock-Tube Problem 

implicit real*8(a-h,o-z) 

dimension u(3,999), un(3,999), ux(3,999), ut(3,999), 

* s(3,999), uxl(3), uxr(3), xx(999) 

data nx/101/, it/100/, ia/l/,dt/.4d-2/,dx/.ld-l/,ga/1.4d()/, 

* ul,u2,u3/.125d0,0.d0,.25d()/,vl,v2,v3/l.d0,().d(),2.5d0/, 

* cr,cl/l.d(),l.d()/ 

c nx must be an odd integer 
c 

nxl = nx 4- 1 
nx2 = nxl/2 
hdt = dt/2.d0 
tt = hdt*dfloat(it) 
qdt = dt/4.d0 
hdx = dx/2.d0 
qdx = dx/4.d0 
dtx = dt/dx 
al = ga - l.dO 
a2 = 3. dO - ga 
a3 = a2/2.d0 
a4 = 1.5d0*al 
do 5 j = l,nx2 
u(l j) = vl 
u(2j) = v2 
u(3j) = v3 
u(l,nx2+j) = ul 
u(2,nx2+j) = u2 
u(3,nx2+j) = u3 
do 5 i = 1,3 
ux(i,j) = O.dO 
ux(i,nx2+j) = O.dO 
5 continue 
c 

open (unit=8,file== , for008’) 
write (8,10) tt,it,ia,nx 
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write (8,20) dt,dx,ga 
write (8,30) ul,u2,u3,cr 
write (8,40) vl,v2,v3,cl 


c 

do 400 i = l,it 
m = nx 4 i - (i/2)*2 
do 100 j = l,m 
w2 = u(2j)/u(l,j) 
w3 = u(3j)/u(l j) 
f 2 1 = -a3*w2**2 
£22 = a2*w2 

f31 = al*w2**3 - ga*w2*w3 
f32 = ga*w3 - a4*w2**2 
f33 = ga*w2 
ut(l j) = -ux(2,j) 

ut(2J) = -(f21*ux(lj) 4- f22*ux(2j) 4- al*ux(3j)) 
ut(3j) = -(f31*ux(l,j) 4- f32*ux(2j) 4- f33*ux(3j)) 
s(lj) = qdx*ux(l j) 4* dtx*(u(2j) 4 qdt*ut(2j)) 
s(2 j) = qdx*ux(2j) 4- dtx*(f21*(u(l j) 4 qdt*ut(l j)) 4- 

* f22*(u(2,j) 4 qdt*ut(2j)) 4 al*(u(3j) 4- qdt*ut(3j))) 
s(3j) = qdx*ux(3j) 4- dtx*(f31*(u(l j) 4 qdt*ut(lj)) 4 

* f32*(u(2,j) 4 qdt*ut(2 j)) 4- f33*(u(3 % j) 4- qdt*ut(3j))) 

100 continue 

if (i.ne.(i/2)*2) goto 150 
do 120 k = 1,3 
ux(k,l)= cl*ux(k,l) 
ux(k,nxl) = cr*ux(k,nx) 
un(k,l) = u(k,l) 
un(k,nxl) = u(k,nx) 

120 continue 

150 jl = 1 - i 4- (i/2)*2 

mm = m - 1 
do 200 j = l,mm 
do 200 k = 1,3 

un(k,j4jl) = 0.5d0*(u(k j) -f u(k,j4l) 4 s(k,j) - s(kj+l)) 
uxl(k) = (un(kj4-jl) - u(k,j) - hdt*ut(k,j))/hdx 
uxr(k) = (u(k,j4l) 4 hdt*ut(kj+l) - un(k,j4jl))/hdx 
ux(kj4jl) = (uxl(k)*(dabs(uxr(k)))**ia 4 uxr(k)*(dabs(uxl(k))) 

* **ia)/((dabs(uxl(k)))**ia 4 (dabs(uxr(k)))**ia 4 l.d-60) 

200 continue 

m = nxl - i 4 (i/2)*2 
do 300 j = l,m 
do 300 k = 1,3 
u(kj) = un(kj) 

300 continue 

400 continue 

c 

m = nxl - it 4 (it/2)*2 
mm = m - 1 
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xx(l) = -0.5d0*dx*dfloat(mm) 
do 500 j = l,mm 
xx(j+l) = xx (j) + dx 
500 continue 

do 600 j = l,m 
X = u(2j)/u(lj) 

y = al*(u(3j) - 0.5d0*x**2*u(l ,j)) 
z = x/dsqrt (ga*y /u ( 1 , j ) ) 
write (8,50) xx(j),u(l j),x,y,z 
600 continue 

c 

close (unit=8) 

10 format(’ t = \f8.4,’ it = ’,i4,’ ia = \i4,’ nx = \i4) 

20 format(’ dt = \f8.4,’ dx = \f8.4,’ gamma = \f8.4) 

30 formatC ul = \f8.4,’ u2 = \f8.4,’ u3 = \f8.4,’ cr = \f8.4) 

40 format (’ vl = \f8.4,’ v2 = \f8.4,’ v3 = \f8.4,’ cl = \f8.4) 

50 formatC x =\f8.4,’ rho =\f8.4,’ v =\f8.4,’ p = \f8.4, 

* ’ M =\f8.4) 

stop 
end 
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Figure I . A surface element on the boundary 
S(V) of an arbitrary space-time volume V. 
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Figure 2. — The SEs and CEs. 
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Figure 3. — Generalized Conservation Elements. 
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Conservation Elements. 



Figure 5. — Unions of Generalized 
Conservation Elements. 
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Figure 8. — GCEs used in the construction of the first, second and 
fourth sets of non-reflecting boundary conditions. 
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Figure 9. — GCEs used in the construction of the third set of 
non-reflecting boundary conditions. 
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(a) Profiles at t=0.2 


Figure 10. = CE/SE solution using the first set of NRBCs 
(Ax = 0.01, CFL= 0.88, a = 1). 
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Pressure Velocity Density 








X 

(c) Profiles at t=0.6 


Fig. 10 (Continued) 
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Pressure Velocity Density 








Figure 11. CE/SE solution using the third set of NRBCs 
with \ = \' =0 (Ax = 0.01, CFL=0.88, a = 1) 
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Figure 12. CE/SE solution at t = 0.4 using the boundary conditions 
Eqs. (5.5) and (5.6) (Ax = 0.01, CFL= 0.88, a = 1) 
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